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ABSTRACT 

We investigate the evolution of low mass (Md/Mb = 0.005) misaligned gaseous discs around 
eccentric supermassive black hole (SMBH) binaries. These are expected to form from ran¬ 
domly oriented accretion events onto a SMBH binary formed in a galaxy merger. When ex¬ 
panding the interaction terms between the binary and a circular ring to quadrupole order and 
averaging over the binary orbit, we expect four non-precessing disc orientations: aligned or 
counter-aligned with the binary, or polar orbits around the binary eccentricity vector with ei¬ 
ther sense of rotation. All other orientations precess around either of these, with the polar 
precession dominating for high eccentricity. These expectations are borne out by smoothed 
particle hydrodynamics simulations of initially misaligned viscous circumbinary discs, result¬ 
ing in the formation of polar rings around highly eccentric binaries in contrast to the co-planar 
discs around circular binaries. Moreover, we observe disc tearing and violent interactions be¬ 
tween differentially precessing rings in the disc significantly disrupting the disc structure and 
causing gas to fall onto the binary with little angular momentum. While accretion from a polar 
disc may not promote SMBH binary coalescence (solving the ‘final-parsec problem’), ejec¬ 
tion of this infalling low-angular momentum material via gravitational slingshot is a possible 
mechanism to reduce the binary separation. Moreover, this process acts on dynamical rather 
than viscous time scales, and so is much faster. 

Key words: accretion, accretion discs - black hole physics - hydrodynamics 


1 INTRODUCTION 

It is widely accepted that most massive galaxies host a supermas¬ 
sive black hole (SMBH) at their centre. Galaxy mergers, expected 
from the hierarchical structure growth scenario based on the A cold 
dark matter (ACDM) cosmological model, then result in the for¬ 
mation of SMBH binaries (Begelman et al. 1980). Both SMBHs 
of such a binary sink towards the galactic centre due to dynamical 
friction and form a hard binary (Merritt 2001). However, most such 
SMBHs appear to be single rather than binary SMBHs, implying 
that SMBH binaries quickly coalesce and merge. One process driv¬ 
ing further binary shrinking are slingshot interactions with stars in 
the ‘loss cone’: those on orbits intersecting with the binary (Saslaw 
et al. 1974). Since the slingshot mechanism ejects the stars, the 
‘loss cone’ needs to be replenished in a relatively short timescale 
in order to shrink the binary all the way down to separations < 10“ 2 
pc where gravitational waves are expected to drive coalescence. For 
spherical collisionally relaxed stellar systems, it is thought that the 
slingshot mechanism stalls well before reaching this separation, re- 

* Email: hal83@le.ac.uk, walter.dehnen@le.ac.uk 
t Einstein Fellow 


suiting in what is known as ‘the final parsec problem’ (Milosavl- 
jevic & Merritt 2003; Berczik et al. 2005). 

Potential stellar dynamical solutions have been sought for gas 
poor systems. Berczik et al. (2006) studied SMBH binaries evolu¬ 
tion in realistic triaxial rotating galaxies and found that the galax¬ 
ies supply stars on centrophilic orbits refilling the loss cone at a 
high enough rate to prevent the SMBH binary from stalling and 
that complete coalescence is achieved in less than 10 Gyr. Khan 
et al. (2013) found that for axisymmetric galaxies with axis ratio 
cl a < 0.8 the hardening rate is 25 times faster than for spherical 
galaxies. Self-consistent A-body simulations of merging galaxies 
containing SMBH found binary hardening rates much higher than 
idealized spherical models and sufficient to shrink the binary to the 
gravitational wave coalescence regime (Khan et al. 2011; Gualan- 
dris & Merritt 2012). 

Interactions with circumbinary gas discs may change the evo¬ 
lution of the SMBH binary. The mass of such a disc is uncertain and 
depends on its formation. A galaxy rich merger can channel large 
amounts of gas towards the centre. If this gas can cool efficiently 
and avoid fragmentation and substantial star formation, a disc with 
mass comparable to that of the binary may form. 

For a prograde disc, spiral density waves in the disc driven by 


2 Hossam Aly, Walter Dehnen, Chris Nixon, & Andrew King 


the outer-Lindblad resonance with the binary transport angular mo¬ 
mentum away from the binary (Goldreich & Tremaine 1979). This 
mechanism is efficient if the disc reaches very close to the binary, 
so that it occupies the resonance, and is sufficiently massive for the 
angular-momentum absorption not to result in an expansion of the 
inner disc edge (Escala et al. 2005; MacFadyen & Milosavljevic 
2008). For a disc with M d = 0.2M b , Cuadra et al. (2009) found that 
binary orbital decay can stall because the disc expands due to ab¬ 
sorption of angular momentum from the binary, severely slowing 
further angular momentum exchange (see also Lodato et al. 2009). 

Apart from the classical density-wave mechanism, the infall 
of gas from the inner edge of the disc into the cavity can be impor¬ 
tant (Roedig et al. 2012; Roedig & Sesana 2014). The binary may 
either eject such infalling gas via a gravitational slingshot whereby 
losing angular momentum and energy, or capture it onto an accre¬ 
tion disc around either component, which adds to the binary angular 
momentum. The binary evolution is determined by the competition 
between these two effects and it remains unclear, which one wins 
in the long term 1 . 

For a retrograde co-planar disc, the lack of orbital resonances 
allows the disc to extend to small radii. This enables the binary 
to accrete or capture material with negative angular momentum 
(Nixon et al. 2011). If M d ~ M b , this may suffice to achieve coa¬ 
lescence (Roedig & Sesana 2014). 

In reality, discs with mass in excess of their aspect ratio times 
the binary mass are gravitationally unstable and hence, due to the 
short cooling time in these discs, fragment and form stars much 
faster than binary coalescence (Gammie 2001; Goodman 2003; 
Levin 2007). The numerical treatment of fragmentation, star for¬ 
mation, and stellar feedback is extremely challenging. In all of the 
aforementioned simulations with such massive discs, these pro¬ 
cesses have simply been suppressed (by assuming slow cooling 
which prevents star formation), overestimating the efficiency of 
disc-driven binary coalescence. Although star formation will rob 
the disc of a significant amount of gas, the newly formed stars 
may still contribute to binary orbital decay (e.g. Sesana et al. 2007, 
2008), though less so than the gas owing to the lack of an efficient 
dissipation mechanism to reduce their pericentres. 

A more likely scenario than binary coalescence driven by the 
interaction with a single massive disc is the repeated interaction 
with low-mass discs resulting from the infall and tidal disruption of 
molecular clouds onto the binary. Such discs are expected to have 
masses 10 5_6 M o typical of molecular clouds, small compared to 
the typical mass 10 6-9 M o of a SMBH binary. Nixon et al. (2011) 
studied retrograde discs of this type and found that they are very ef¬ 
ficient in reducing the binary angular momentum through accretion 
of gas with negative angular momentum onto the secondary black 
hole. This enhanced accretion onto the secondary black hole in¬ 
creases the binary’s eccentricity, decreasing the pericentre distance 
in the process, and coalescence is achieved when a mass compara¬ 
ble to the secondary black hole has been accreted. 

If accretion events in galactic nuclei are chaotic and randomly 
oriented (King & Pringle 2006, 2007; King et al. 2008), we ex- 

1 This effect was present in the simulations of Cuadra et al. (2009) but did 
not effectuate significant binary evolution. On the other hand, based on an 
extrapolation to 50 times longer than actually modelled, Roedig & Sesana 
(2014) claim efficient binary shrinking. However, since the infall of gas de¬ 
pends on the disc structure at its inner edge, this result is very sensitive to the 
thermodynamical treatment. Roedig et al. (2012), for example, found that 
for isothermal instead of adiabatic gas with an imposed standard ,6 cooling 
prescription, the binary orbital decay can be significantly reduced. 


pect the formation of misaligned circumbinary discs around SMBH 
binaries. In the case of a circular SMBH binary, the interaction 
between the misaligned disc and the binary is similar to Lense- 
Thirring precession on an accretion disc around a spinning black 
hole (Bardeen & Petterson 1975; Pringle 1992; Scheuer & Feiler 
1996). King et al. (2005) showed that the induced differential pre¬ 
cession will cause a misaligned disc to counter-align with the black 
hole spin provided 


where J d and / h are the disc and black hole angular momenta, re¬ 
spectively; and 6 is the angle between them. The disc will co-align 
with the black hole spin if this relation is not satisfied. Nixon et al. 
(2011) showed that the same analysis applies to the case of a mis¬ 
aligned disc around a binary (though the precession rate is slightly 
different). Thus, if counter-alignment is stable (Nixon 2012), this 
mechanism can provide a solution to the final parsec problem by 
supplying retrograde discs to achieve coalescence. 

Recently, Nixon et al. (2013) performed 3-D hydrodynamical 
simulations of circumbinary discs around a circular binary for var¬ 
ious tilt angles 6. In addition to co- and counter-alignment, they 
found that in many cases the discs is torn into distinct rings which 
precess almost independently (Nixon et al. 2012). The precessing 
rings, which have partially opposed angular momentum, may in¬ 
teract causing partial cancellation of their angular momenta and 
thus gas infall close to the binary. This disc tearing significantly 
increases the accretion rate and may play an important role in pro¬ 
moting the binary final coalescence. 

Those studies considered the case of a circular binary interact¬ 
ing with a circumbinary disc, when disc precession is only around 
the pole of the binary plane. In this study, we consider the more 
general situation of an eccentric binary. For a SMBH binary formed 
via a galaxy merger, we expect high eccentricities in many cases 
(Aarseth 2003; Khan et al. 2011; Wang et al. 2014). Moreover, ret¬ 
rograde accretion onto a circular binary naturally results in eccen¬ 
tricity growth as discussed earlier. One important effect of binary 
eccentricity is to make the time averaged binary potential triaxial 
rather than axisymmetric as for a circular binary. Previous studies 
have shown that misaligned discs in triaxial galaxies can precess 
around both the major and the minor axes (Steiman-Cameron & 
Durisen 1984; Thomas et al. 1994). 

The paper is organised as follows. In Section 2 we present an¬ 
alytic results for a simple orbit-averaged model for the binary-disc 
interaction up to quadrupole order. Section 3 describes the setup 
of our 3D hydrodynamical simulations, the of which results are 
presented in Section 4 and discussed in Section 5. Finally, we sum¬ 
marise and conclude in Section 6. 


2 BINARY-DISC QUADRUPOLE INTERACTION 

The dynamics of a circumbinary gaseous ring orbiting an eccentric 
binary is not analytically treatable, even without considering any 
dissipation. However, useful insight can be obtained by (i) truncat¬ 
ing the binary gravitational potential at quadrupole order, (ii) as¬ 
suming that the ring is circular, (iii) time-averaging over the binary 
orbit, and (iv) neglecting dissipation. Assumptions (i) and (ii) are 
valid as long as the ring is sufficiently distant from the binary, while 
assumption (iii) requires that orbital resonances between ring and 
binary are not important. 

The monopole of the gravitational interaction results in Kep- 



lerian motion of the ring around the binary centre of mass, while 
the quadrupole describes the lowest-order deviation of the binary 
from a central point mass. 

Recently, Naoz et al. (2013) have used Hamiltonian perturba¬ 
tion theory to obtain the equations for the secular evolution of a 
hierarchical triple up to octopole order. For a circular outer binary, 
their results are equivalent to the situation of a circular circum- 
binary ring. We now summarise the relevant relations (obtained in 
Appendix A with Newtonian dynamics, but otherwise equivalent to 
those of Naoz et al.) in terms of vectors rather than orbital elements 
to describe the system. 

The binary is parametrised by its mass ratio q = m 2 lm\ < 1, to¬ 
tal mass M-m\ + m 2 , semi-major axis a, specific angular momen¬ 
tum h, and eccentricity vector e. Let R = x\-x 2 the instantaneous 
binary separation vector, then 


, m\ m 2 q „ • 

h- —X 1 XX 1 + — x 2 x x 2 —- —RxR and 

M M (1 +q) 2 

Rx(RxR) h 

6 = —GM -* 


( 2 ) 

(3) 


The vector e is conserved for the binary orbit and points from the 
centre of mass to peri-apse (hence is always orthogonal to h). Its 
magnitude is the orbital eccentricity and is related to that of h by 


h 2 (l+q) 4 = q 2 GMa(l-e 2 ). 


(4) 


2.1 Ring evolution 

The circular circumbinary ring is parametrised by its mass m, ra¬ 
dius r, and pole /. The latter is the unit vector in direction of 
the ring’s specific angular momentum /, which has amplitude / = 
^G(M + m)r. The ring radius must satisfy r > a(\ + e)/(l + q) for 
the quadrupole-approximation to be valid (and in order to avoid 
collision with a binary component). Note that the tilt angle 6 of the 
ring with respect to the binary satisfies cos 6 - l-h. 

The quadrupole interaction energy between binary and ring, 
averaged over both the binary orbit and the ring, is 

mco 2 a 2 cj ro 0 a a 0 -i 

(E br ) = - g * [6e 2 -1 - I5e 2 (l■ e) 2 + 3( 1 - e 2 )(l■ h) 2 ] (5) 

with to = \/G( M + )h)/ r‘ the orbital frequency of the ring, in agree- 
ment with equation (22) of Naoz et al. (2013). The time-averaged 
binary quadrupole torques the ring according to 

/ = 0x/ (6) 

with the vector 

0= 4iTT^$[ 5e2(i -^“ (1 “ e2)(24) ^- (7) 

From equation (6) we have / • / = 0, i.e. / = 0 = r and the ring is 
merely precessing (this is no longer true at octopole and higher 
order, when / ± 0, see Naoz et al. 2013). 

Since tumr 2 0 = d(E br )/dl, equation (6) implies 

0. (8) 

CIf dl 

Thus, (in the assumed approximation) no energy is exchanged be¬ 
tween binary and a single ring, but only angular momentum (if the 
binary interacts with several rings, the individual interaction ener¬ 
gies with each ring are no longer conserved). 
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2.2 Binary evolution 


The torque of the binary from the ring can be worked out analo¬ 
gously to that of the ring from the binary. After averaging over the 
binary orbit, we obtain 


w 

h = - 0x/. (9) 

M 

In particular, the total angular momentum, Mh+ml , is conserved at 
quadrupole order. For the case w«M considered here, the orien¬ 
tation h only varies slightly even if the disc orientation / undergoes 
large changes. 

For a circular binary 0 is parallel to h such that hh- 0, i.e. 
h = \h\ is conserved and the binary is merely precessing (with an 
amplitude that is smaller than that of the disc by a factor m/M). 
This fact together with conservation of total angular momentum 
was the basis of the analysis by Nixon et al. (2011). 

For an eccentric binary, the evolution of h is not simply a pre¬ 
cession and h not conserved. Instead, we find 


15 a) 2 m e 2 h 
4 O-M y ] _ e i 


Cl-e) (Ik) 


( 10 ) 


with f2 = -y GM/a 3 the binary orbital frequency. Thus, h remains 
unchanged only if e = 0 (circular binary), or if / is perpendicular 
to either e or k , i.e. if either e or k are in the ring plane. Other¬ 
wise, h oscillates, since l • k oscillates around zero under the ring 
precession. 

The change of the eccentricity vector is 

4 = iS¥ eVrr ^{[ 2_(2 * )2_5(i ' g)2 ] k+ 

+ (l-k)(l-h)h + 5(l-e)(l-k)e} (11) 

and the corresponding change in eccentricity 

e = L ( ddd e Vl — <? 2 (/ ■ e )(/ • k). (12) 

in agreement with equation A34 of Naoz et al. 2013, but also with 
equation (10) in conjunction with equation (4). In addition to the 
precession of the orbital plane and the oscillation of the eccentricity 
(both already described by equation 10), the binary also undergoes 
apsidal precession with rate 

Vi-e 2 \l-(l hf-5(l-e) 2 ], (13) 

which is prograde for near-planar disc orientations (when \l-h\ ~ 1), 
but retrograde for near-polar discs (when \l-e\~ 1). 


2.3 Ring precession 

The rates of change of the directions of the binary and ring angular 
momenta satisfy 


ml 

Mh 


(14) 


and a similar relation holds for \de/dt\. Thus, as long as m « M, 
the binary orientation changes only very little and/or much more 
slowly than that of the ring (except for extreme binary eccentric¬ 
ities when Mh oc ^1-e 2 can be small). We therefore consider in 
this subsection the limit m/M ^ 0 when the binary orientation and 
eccentricity are constant. 
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Figure 1. Precession paths for the direction / of the angular momentum of a dissipation-less circular ring of negligible mass orbiting a binary with eccentricity 
e as indicated. The binary orbits counter clockwise in the plane perpendicular to its specific angular momentum h with peri-apse in the direction e. For a 
circular binary (e = 0, left), / always precesses around h in a retrograde sense. For eccentric binaries, prograde polar precession (blue) around e, the long axis 
of the time-averaged binary potential, is also possible. The regions of polar and azimuthal precession are separated by two great circles (black). The four ring 
orientations / = ±h and l-±e are stable (non-precessing), while the orientations l = ±k are unstable. Dissipation would damp the precession and eventually 
align the ring with one of the four stable orientations. In case of a massive ring, the binary orbit evolves too: the vectors h and e oscillate and precess, and e 
and k rotate around h. 


Then, equation (8) implies that the ring precesses along curves 
of constant (E hr ). Isolated minima and maxima of (E br ) denote sta¬ 
ble, non-precessing ring orientations. In the presence of dissipa¬ 
tion (due to viscosity in the disc), these orientations are attractors, 
i.e. the dissipative damping of the precession eventually aligns the 
pole / of the ring with the extrema of (E br ) (Steiman-Cameron & 
Durisen 1984). For any e < 1, the orientations l = ±h are isolated 
minima of (E b r ) and correspond to co-planar ring orientations ei¬ 
ther co- or counter-rotating with the binary. 

For a circular binary (e = 0), these are the only stable orien¬ 
tations, but all polar orbits (6 = 90°) maximise (E br ). © is parallel 
to h and ring precession is circular: / describes a circle around ei¬ 
ther of the stable orientations, see also the left plot in Fig. 1. The 
precession rate is lower than the orbital frequency by the factor 
3qa 2 cos6/4r 2 (l +q) 2 . This is the situation previously studied by 
Nixon et al. (2011). We now turn to the more general case of an 
eccentric binary. 

For e > 0, the orientations l-±e are maxima of (E br ), corre¬ 
sponding to polar rings (with opposite senses of rotation) around 
e. For e < 1, l—±k are saddle points and correspond to polar rings 
around 

k = hxe , (15) 


the intermediate axis of the time-averaged binary potential. 
These latter ring orientations are unstable, i.e. small deviations will 
result in precession around either of the four stable orientations. 
For 0 < e < 1, ring precession is never circular: / describes a curve 
elongated towards the unstable orientations, rather than a circle. 
Azimuthal and polar precessions are retrograde and prograde, re¬ 
spectively. See Fig.l for a visualisation of the precession paths. 

The regions of polar and azimuthal precession are separated 
by the contours of (E br ) passing through the saddle points. These 
separatrices are circular and shown as black in Fig. 1. The fraction 
of ring orientations undergoing polar precession is 


- cos 

71 


\-6e 2 

l+4e 2 ' 


(16) 


At small e, this grows linearly (oc V20 e/n) with eccentricity (see 
Fig. 2). Azimuthal and polar precession are equally likely for e = 
6“ 1/2 w 0.408. 



Figure 2. Percentage of ring orientations undergoing polar precession as a 
function of binary eccentricity. 


3 SIMULATION SETUP 

We perform a set of 3-D Smoothed Particle Hydrodynamics (SPH) 
(Gingold & Monaghan 1977; Lucy 1977) simulations of geometri¬ 
cally thin accretion discs with different initial misalignment around 
an eccentric binary. We use a range of different binary eccentrici¬ 
ties e - 0, 0.3, 0.6, and 0.9. The disc setup is very similar to that 
used by Nixon et al. (2013): the disc is initially flat and extends 
from an inner radius of 2 a to an outer radius of 8a with an inner 
thickness H/R = 0.01. We use a disc viscosity coefficient (Shakura 
& Sunyaev 1973) a = 0.1 which we setup using an appropriate SPH 
artificial viscosity coefficient a A v corresponding to our resolution 
(Lodato & Price 2010). All simulations start with 4 million SPH 
particles, while the the binary is modelled using two equal mass 
sink particles with accretion radius of 0.05a. The disc initial surface 
density follows the profile 'LocR~ 3/2 , and we use a locally isother¬ 
mal equation of state with sound speed c s ocR~ 3/4 . These choices 
ensure a uniform vertical resolution (and hence uniform physical 
viscosity, see Lodato & Pringle 2007). We assume a disc mass of 
M d /M b = 0.005 < H/R, which ensures that disc self-gravity is not 
important (we do not include gas self-gravity in our simulations, but 
we do self-consistently include the back-reaction from the gas on 
the binary). The simulations were performed using our own code 
(Dehnen & Aly 2012), which implements an SPH scheme very 
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6= 30° 


0 = 45° 


<9 = 60° 


0 = 80° 


0 = 90° 


e = 0 


e = 0.3 


e = 0.6 


e = 0.9 



Figure 3. Density rendering of the simulations after t = 600 (« 95 binary orbits) for different binary eccentricities and initial misalignment angles as indicated. 
The projections are along the intermediate binary axis k with the angular- momentum vector h pointing upwards and the eccentricity vector e to the right. 


similar to that used by Nixon et al. (2013), and we verified that 
our results for e - 0 agree with theirs. The disc has initial angular 
momentum direction 


/= sin#cos0£ + sin#sin0& + eos#/*, (17) 


where 0 is the twist angle of the disc. We ran a total of 118 simula¬ 
tions for e = 0, 0.3, 0.6, and 0.9; 6 = 0°, 10°, 30°, 45°, 60°, 80°, 90°, 
100°, 120°, 135°, 150°, 170°, and 180°; 0 = 0°, -45°, and -90.°. 
In the next section 0 will be taken to be zero whenever it is not 
specified, we discuss the effects of varying 0 separately. 

We point out that the choice of a locally isothermal equation of 
state implies that the disc instantly radiates away all the heat gained 
from the viscous dissiaption and shocks. This is justified if the the 
cooling time is much shorter than the precession time. When this 
assumption does not hold, the disc thickness will increase and will 
be more able to resist breaking. We leave more advanced thermo¬ 
dynamic treatment to future investigation. 


4 SIMULATION RESULTS 

Fig. 3 shows snapshots after « 95 binary orbits for the twenty sim¬ 
ulations with initial tilt angles 6 = 30°, 45°, 60°, 80°, and 90° ini¬ 
tial binary eccentricities e - 0, 0.3, 0.6, and 0.9. As expected from 
the analysis in Section 2, the disc precesses around h for circular 
and low-eccentricity binaries, while for very high eccentricities the 
precession is predominantly around e. In almost all cases, the disc 
breaks into distinct rings, which interact with each other and, de¬ 
pending on the details of each case, result in either co-, counter-, or 
polar-alignment of the disc. In some cases the interaction between 
independently precessing such rings is very violent and disruptive, 
leading to ejection of gas. We now visit each possible outcome in 
detail. 

4.1 Polar alignment 

Nixon et al. (2011) showed that for a circular binary, where the bi¬ 
nary induced precession is only around h , the disc eventually co- or 
counter-aligns with the binary orbital plane depending on the disc 
angular momentum and its initial misalignment angle. Our analy¬ 
sis in Section 2 suggests that for an eccentric binary the precession 
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R/a 



R/a 

Figure 4. Evolution of tilt profiles for e = 0 (top) and e - 0.9 (bottom) discs 
with initial tilt 6 = 60° at t = 0, 100, 200, 300, 400, 500, and 600 (see legend) 
in code units. 


will be around h or e. In the latter case, dissipation results in polar- 
alignment. 

Fig. 4 shows the time evolution of the tilt profiles 6{R) for two 
simulations with initial 6 = 60° but either e = 0 (top) or e - 0.9 (bot¬ 
tom). For the circular binary case, the inner part of the disc even¬ 
tually co-aligns with the binary (6 —» 0°), as expected. In contrast, 
for the highly eccentric binary we see the disc aligning in a polar 
configuration with respect to the binary angular momentum vector 
(6 —» 90°). The variations in 0 as function of both time and radius 
are caused by the precession around e. 

This is more evident from Fig. 5, where we plot the orientation 
/ of the angular momentum in annuli of the disc at t = 100 for four 
simulations with different binary eccentricity but identical initial 
disc orientation at 6- 60°. We see that the discs in our simulations 
closely follow the predicted precession paths especially in the outer 
parts of the disc. The inner parts of the disc, which have higher 
precession rates, dissipate faster and start to align with the h or e as 
expected. 


4.2 Violent ring interactions 

Our simulations starting from discs misaligned to both the h and e 
show rather violent gas dynamics. The radially differential binary 
torque tears the disc and causes the formation of separate rings. 
These rings are mutually misaligned and start to interact with each 
other, presumably because they gained some eccentricity from in¬ 
teractions with the binary. The ring interactions cause partial can¬ 
cellation of angular momentum and hence a significant increase in 


the accretion rate. This is identical to the picture reported by Nixon 
et al. (2013) for circumbinary discs around circular binaries. 

However, for very high eccentricities we find disc tearing to be 
much more violent and lead to a different evolution from that for 
circular and low-eccentricity binaries. There are two reasons for 
this difference: first the precession rate increases with eccentricity; 
second, the low-angular-momentum gas resulting from the inter¬ 
actions and falling onto the binary will align to polar orientation 
rather than a prograde or retrograde orientation as in the case of a 
near-circular binary. This allows this highly eccentric low-angular- 
momentum gas to come very close to the binary without suffering 
a lot of accretion. This non-circular gas in the central zone interacts 
with the outer disc further increasing its orbital eccentricity, throw¬ 
ing more gas to the centre, and promoting more interaction. This 
run away effect is shown in Fig. 6 and can also be seen in the bot¬ 
tom left panel of Fig. 3. Eventually, this process sends an increasing 
amount of gas plunging onto the binary on almost radial orbits that 
can reach very close to the binary, avoiding significant accretion, 
and receiving energy kicks from one of the binary components in 
a manner very similar to the slingshot mechanism. Some of this 
gas will get ejected producing outward, almost radial, streams that 
can act as a possible observational signature of a highly eccentric 
SMBH binary. 

Fig. 6 shows density rendering (left panels) and particle plots 
coloured by eccentricity magnitude (right panels) of 5 snapshots 
for the simulation with e = 0.9 and initially 6 = 150° at times t = 
0, 200, 400, 600, and 800. We can see that the amount of chaotic 
gas resulting from the ring interactions keeps increasing during the 
simulation. The outward streams of gas resulting from the slingshot 
mechanism are very clear. 


4.3 Precession rate 

In order to provide a quantitative comparison between the predic¬ 
tions of our analytical model in Section 2 and the results obtained 
from the simulations, we plot in Fig. 7 the analytical precession rate 
0 derived from equation (7) for a disc with an initial misalignment 
of 6- 60° around binaries with four different eccentricities along 
with the equivalent precession computed from the simulation and 
averaged over 10 binary orbits starting from t = 50. We find that the 
simulations agree quite well with the predicted precession rate at 
radii > 2.5 R/a. For discs around eccentric binaries, we observe os¬ 
cillations on binary orbital timescales and a good agreement is only 
found when the precession rate is averaged over a few binary orbits. 
We note that only a modest agreement is to be expected since our 
model ignores dissipative effects, contributions from higher than 
quadrupole order, and orbital resonances. 

4.4 Non-zero initial disc twist angle 

So far, all the results shown here are for twist angles 0 = 0°, i.e. ini¬ 
tially the disc line of nodes with the binary plane is the k direction: 
/ is tilted towards e. In general, however, we should expect any disc 
orientation, i.e. non-zero 0. 

In Fig. 8 we present snapshots for simulations with binary ec¬ 
centricity e - 0.9, initial twist angles 0 = 0°, 45°, and 90°, and ini¬ 
tial tilt angles 6 = 30°, 45°, 60°, 80°, and 90°. For 0 = 90°, we only 
observe azimuthal precession, akin to the circular binary induced 
precession. This confirms our prediction since for that case / and e 
are always orthogonal, causing the first term in equation (7) to van¬ 
ish, and we are left with only azimuthal precession. For 0 = 45°, 
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Figure 5. Projections of angular momenta of the radially binned disc in our simulations compared to the analytical precession paths shown in Fig. 1. Solid 
circles represent seven radial bins of the disc ranging from R = 1 (light gray) to R = 8 (dark gray) at t- 100 for simulations of different eccentricity (as 
indicated) but the same initial tilt 6 = 60°. Obviously, the inner disc precesses faster, nicely following the theoretical precession paths. The innermost parts of 
the discs around eccentric binaries start to align with the stable polar orientation. 


we find the same trend discussed earlier, i.e. experiencing either 
polar or azimuthal precession, or violent ring interaction. In Fig. 9 
we compare the precession paths of all three 0 values for the case 
of e = 0.9 and 6 = 60° to the analytical contours. Similar to Fig. 5, 
we see the simulations closely follow the analytical contours apart 
from the innermost parts where disc disc breaking and aligment are 
dominant. This strongly suggests that our analysis above for 0 = 0° 
carries over to the general case. 


5 IMPLICATIONS FOR THE FINAL-PARSEC PROBLEM 

The solution suggested by Nixon et al. (2011) to the final parsec 
problem requires the binary to accrete negative angular momen¬ 
tum from a retrograde disc, which gradually increases the binary 
eccentricity until coalescence is achieved via energy losses to grav¬ 
itational radiation at pericentre. Nixon et al. (2011) showed that for 
a circular binary counter-alignment of randomly oriented accretion 
events can provide a continuous supply of the required retrograde 
discs. In particular, they showed that for cases where J b > 2/ d , all 
accretion events with initial misalignment of# > 90° will result in 
a retrograde disc, i.e. roughly half of randomly oriented accretion 
events will counter-align with the binary as long as the binary dom¬ 
inates the angular momentum of the system. 

Our results somewhat change this picture. As the binary ec¬ 
centricity increases (due to retrograde accretion as in Nixon et al. 
or earlier stellar dynamical processes) disc counter- (and co-) align¬ 
ment becomes ever less likely at the expense of polar alignment. 
The subsequent accretion of such polar discs merely rotates the an¬ 
gular momentum vector of the binary presumably hardly affecting 
the binary eccentricity. Thus simply retrograde gas accretion ap¬ 
pears less viable a solution to the final parse problem. 

There are however, still several ways gas can solve this prob¬ 
lem. First, a single massive retrograde accretion event may, in prin¬ 
ciple, supply enough negative angular momentum to complete the 
binary merger. However, for a massive disc self-gravity becomes 
important, likely causing clumping and star formation, which re¬ 
duces the amount of gas that can be accreted. Moreover, a single 
massive retrograde accretion event may be not be sufficiently likely 
to explain the coalescence of all SMBH binaries (which form with 
each major merger of massive galaxies). 

A more intriguing possibility involves more violent gas dy¬ 
namics. We showed that, in many cases, the disc does not smoothly 
align, instead the strong differential precession (in particular for 
misaligned discs around eccentric binaries) leads to tearing of the 


disc into separate mutually misaligned rings. In the inner disc close 
to the binary, the gravity of the binary cause these rings to become 
eccentric such that they inevitably interact with each other and with 
the outer disc. These interactions cause further eccentricity growth 
on a dynamical time scale and eventually result in plunging gas 
infall. Some of this infalling gas will be accreted by either binary 
component. This will change the binary angular momentum, but 
may not reduce its absolute value, depending on the orientation and 
in contrast to the situation with pre-dominantly retrograde accretion 
(Nixon et al. 2013). 

If the infalling gas evades this fate, it will most likely get 
ejected from the binary via a slingshot interaction. This in turn re¬ 
duces the binary separation in much the same way as the ejection of 
penetrating stars, thus exactly as required to solve the final parsec 
problem. Indeed, we find in our simulations which undergo violent 
gas dynamics not only significant gas accretion but also a shrinking 
of the binary orbit. 

Clearly, this violent interaction and accretion processes are 
rather complex and chaotic and certainly not well resolved or ad¬ 
equately modelled in our simulations. Nonetheless, what our sim¬ 
ulations quite clearly show is that such violent gas-dynamical pro¬ 
cesses are inevitable if the gas is initially misaligned with the bi¬ 
nary, in particular if the binary is eccentric. We leave a more de¬ 
tailed investigation into the binary orbital evolution in this chaotic 
environment for a future study. 


6 CONCLUSIONS 

We have studied the interaction of an eccentric binary with a 
gaseous disc initially misaligned with the binary angular momen¬ 
tum. Such a configuration should occur naturally from the infall 
and subsequent circularisation of gas into the inner few parsec of 
a merger remnant still hosting a supermassive black hole (SMBH) 
binary (e.g. Dunhill et al. 2014). The binary exerts a torque on the 
disc, resulting in disc precession and, due to viscous dissipation, in 
eventual alignment of the disc with the binary. In case of a circular 
binary, this alignment is always co-planar, resulting either in a pro- 
or retro-gradely rotating circumbinary disc (Nixon et al. 2011). 

We find that in the general case of an eccentric binary po¬ 
lar alignment also occurs, when disc angular momentum is aligned 
with the binary peri- or apo-apse direction. The binary torque on the 
disc can be quite accurately understood analytically assuming an 
orbit-averaged binary potential to quadrupole order (see Section 2). 
The fraction of initial disc orientations which give rise to polar 
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Figure 8. Density rendering of the <? = 0.9 simulations at t = 600 (« 95 orbits of the binary) projected on the x-z plane with different values for 0 and 6 as 
indicated in the figure. 


alignment grows with binary eccentricity, reaching 0.5 at e « 0.4. 
The precession paths (neglecting dissipation) are not circular, but 
elongated towards the intermediate axis of the orbit-averaged bi¬ 
nary. 

The prospect of accretion onto the binary from a polar instead 
of co-planar disc impedes the solution (proposed by Nixon et al. 
2011) to the final-parsec problem for coalescing a supermassive 
black-hole binary. In that picture consecutive randomly oriented 
accretion events lead to the formation of either pro- or retro-grade 
co-planar circum-binary discs. While accretion from the former is 
largely suppressed (by orbital resonances as discussed in the intro¬ 
duction), accretion from the latter reduces the binary angular mo¬ 
mentum and drives it to larger eccentricities. However, at large ec¬ 
centricities polar disc orientations dominate, when accretion (not 
resolved in our simulations) has presumably little effect on the bi¬ 
nary orbit (since the accreted angular momentum is perpendicular 
to that of the binary). Thus, eccentricity growth via accretion is 
likely to be significantly reduced well before gravitational wave 
emission can take over as driver for coalescence. 

However, in many of our simulations, in particular for larger 
binary eccentricity and stronger initial misalignment, the disc does 
not smoothly align, but is torn into separate mutually misaligned 
rings. This process was already reported by Nixon et al. (2013) for 
a circular binary and can be understood by the radially differential 
binary torque, which overcomes the adhesive effect of gas viscos¬ 
ity. The prominence of tearing with binary eccentricity and initial 
disc misalignment can be understood as consequence of the larger 
binary torque in these cases. 

The subsequent evolution of these gas rings can be rather 
chaotic and is not quite adequately modelled in our simulations. 


However, some basic results appear to be robust. The innermost 
rings are sufficiently perturbed by the binary to acquire some or¬ 
bital eccentricity. This in turn inevitably leads to interactions be¬ 
tween the rings, resulting in partial cancellation of their angular 
momenta. This process is more prominent in more eccentric bi¬ 
naries, because the stronger binary torque results in larger mutual 
misalignment between adjacent rings. The cancellation of angular 
momentum of the rings will increase their eccentricity, providing 
a positive feedback loop and hence a run-away process, eventually 
resulting in gas plunging onto the central binary. This material may 
be accreted onto either hole, but when coming from a near-polar 
orientation, this will hardly help with the final-parsec problem, as 
explained above. 


Alternatively, the infalling gas, which for a highly eccentric 
binary can come much closer to the binary whilst avoiding accre¬ 
tion, may get ejected from the binary via a gravitational slingshot 
interaction with one of its components. This also helps to solve the 
final-parsec problem, though this time by reducing its semi-major 
axis. This is similar to the stellar-dynamical process of shrinking 
the binary orbit via ejection of stars penetrating into the binary. 
The difference is that the total amount of stars in the ‘loss cone’ 
(whose orbit carries them into inner parsec) is limited and cannot 
be easily re-filled, while gas being dissipative and collisional by 
nature may provide a better agent. This is particularly the case at 
the parsec scale where the SMBH dominates the dynamics and by 
its gravitational torques shepherds some gas into the loss cone. 
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Figure 6. Density rendering (left panels) and particle plots coloured by ec¬ 
centricity magnitude (right panels) of 5 snapshots of the e = 0.9 6 = 150° 
run at times (from top to bottom) t- 0, 200, 400, 600, and 800. 



b=_i_i_i_I_i_i_i_I_i_ ■ ■ -H 

2 4 6 8 
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Figure 7. Comparison between the disc precession rate measured from the 
simulation (solid) with initial 6 = 60° and e = 0, 0.3, 0.6, and 0.9 and our an¬ 
alytical model (dotted) of equation (7). For each case we plot the dominant 
component of 0, i.e. &h for e = 0, 0.3 and & e for e = 0.6, 0.9. 
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APPENDIX A: BINARY-DISC QUADRUPOLE 
INTERACTION 

Here, we give the details of the analysis leading to the results re¬ 
ported in Section 2. Our results, obtained via Newtonian dynamics, 
agree with the more general results of Naoz et al. (2013), obtained 
via Hamiltonian perturbation theory, for a circular ring and ignor¬ 
ing octopole terms. 

The three unit vectors h, e, and k-hxe are conserved under 
the binary motion and are mutually orthogonal, such that 

hihj + efij + kikj = 5ij . (A 1) 

The binary components are at positions 

X! = -^-R, x 2 = --^—R (A2) 

\+ q 1 + q 

with 

R = a(cosq-e)e + a Vl-e 2 sin 77 A:. (A3) 

Here, q is the eccentric anomaly, which is related to the mean 
anomaly £ via 

^ = 77 - e sin 77 , (A4) 

such that d£ = (1 - e cos 77^77 and an orbit average becomes (•) = 
(2^)- 1 C-d -ecos 77 ^ 77 . When orbit-averaging RiRj, the cross term 
between e and k averages to zero and 

(RiRj)= ±tf 2 [(l+4e 2 )e ; e; + (l-e 2 )£,fc,] (A5) 

= ja 2 [Se 2 eiej + (l -e 2 )(5,y —, (A 6 ) 

where the second form follows from eliminating k t kj in favour of 
hihj with the help of the identity (Al). From this result, we can 
work out the orbit-averaged trace-free specific quadrupole moment 


of the binary as 

Q y = M~ l \mi{x u xy-\x\6 i j) + m 2 {x 2i x 2 j-\x\6 i j)\ (A7) 

= [{RiRj) - ^RkRASij] (A8) 

= Q7^[(F e2 F + l deiej-^-e^khj]. (A9) 

We will also need the orbit average 

( R t RjR k ) - ^a 3 Cle V l-e 2 (eikje k + k t eje k - 2^4) • (A 10) 

with Q = 4 GM/a 3 the binary orbital frequency. 
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A1 Ring evolution 


A ring particle at position r experiences the orbit-averaged 
quadrupole potential of the binary 


^ ^ 3GM r, 

(O b )(r) = -^-r-Q-r. 


(All) 


Averaging over the ring, we obtain 

{r i r j ) = \r\8 i j-U j ), (A12) 

such that the trace-free specific quadrupole moment of the ring is 


°l ij~ r \ 6 


(A13) 


The quadrupole interaction energy between binary and ring, aver¬ 
aged over the the binary orbit and the ring, is then 

3 GMm 


<£br> 


2 r 5 


-tr(e) 


(A 14) 


with interaction tensor O = q Q, which has components 

[(e -<? 2 )<5,; + f - 5(1 -e 2 )hihj ~ (A15) 

(± -3e 2 )l i ! j -^e 2 (i-e)! i e j + §(1 -e 2 )(l■ h)hhj\. 


a 2 r 2 q 

Qii ~ 6(1 +q) 2 


Taking its trace in equation (A 14) we obtain equation (5). The orbit- 
averaged torque on the ring also involves the interaction tensor. Us¬ 
ing index notation, we have 

3GM -e ijk Ojk■ (A16) 


/ / 8 

Note that only the anti-symmetric part of O contributes to the 
torque. Inserting equation (A 15), we find that we can write this 
as / = 0 x / with the vector 

3 cjq a 2 ■ 


0 = 


■[5e 2 (l-e)e-{\-e 2 )(l-H)h], 


4(1 +qf r 2 

where a> = + m)/t s is the ring angular frequency. 


(A17) 


A2 Binary evolution 


The torque of the binary from the ring can be worked out analo¬ 
gously to that of the ring from the binary. The quadrupole potential 
due to the ring at the binary is 

3 Gm 

®r = —(A18) 

Adding the torque from each binary component and averaging over 
the orbit, we find 
3 Gm 

hi = —^s ijk Q kj . (A 19) 

In particular, the total angular momentum, Mh + m/, is conserved 
at quadrupole order. Together with the precession of the ring, this 
implies that the evolution of h is not simply a precession and that 
h = \h\ is in general not conserved. Instead, we find 


o/m 15 e 2 h 
Q.M 4 yi _ e 2 




(A20) 


Thus, h remains unchanged only for a circular binary, h = 0 for e = 0 
or if l is perpendicular to either e or k. Otherwise, h oscillates, since 
/ • k oscillates around zero under the ring precession. 

The change of the eccentricity vector is 


2 R(R R) - R(R R) - R(R R) 

GM 


(A21) 


with 

dO r _dO r _3 Gm 
dx 2 dx\ r 5 


(A22) 


Inserting (A22) into (A21) and orbit averaging (using equa¬ 
tion A 10 ), we find 


e - 3- e 

Q.M 


U [ke q e-ek c\ e-\k c\] 


= 3 ~e sTTe 2 [({-d-e) 2 )k + (le)d-k)e+{d- 
= J^ e ^ I {[ 2 -( i -h) 2 -5(le) 2 }k + 

+ (l-k)(l-h)h + 5(l-e)(l-k)e}. 

L -e^ll -e 2 (l• e)(l• k) 


(A23) 

(A24) 

(A25) 


From this, we obtain 
15 (j) 2 m 

in agreement with equation A34 of Naoz et al. (2013), but also with 
equation (A20) and (1 -e 2 )h = -hee (from equation 4). 


(A26) 


A3 Ring precession 

If the mass of the ring is negligible compared to that of the binary, 
we can approximate the binary orientation as fixed and the vectors 
h , e , and k as constants. In this case, the evolution of the ring orien¬ 
tation allows some further analytical treatment. 

Since 0 is parallel to d(E hr )/dl , precession is along lines of 
constant (E br ). This gives the equation 

C = (1 - e 2 )(l • h ) 2 - 5 e 2 (l • ef (A27) 


with constant C for the precession paths. C = 0 corresponds to the 
contour of (E br ) through the unstable orientations l = ±k. Hence, 
this contour separates the regions of polar and azimuthal preces¬ 
sion. The r.h.s. of equation (A27) can be written as (/ • u{) (1 • u 2 ) 
with 

U U2 = ± s[5e. (A28) 


Thus, the separatrices are great circles with poles u 1>2 . The fraction 
of ring orientations undergoing polar precession is 


- cos - u 2 )= — cos' 

n n 


\-6e 2 
1 +4e 2 


(A29) 


At small e , this grows linearly (oc V20^/7r) with eccentricity. Az¬ 
imuthal and polar precession are equally likely for u\ u 2 = 0 , which 
occurs at e - 6~ l/2 « 0.408. 

If the constant C in equation (A27) is positive (negative), we 
have azimuthal (polar) precession. This equation for / has the para¬ 
metric solutions for the precession paths (see also Fig. 1) 


/ -8 = 


-C 


1 +4e 2 


cos£, Ik- 




-e 2 -C 


1 - 


sin^ 


for 0 < C < 1 -e 2 (azimuthal precession), and 

1 f / 5e 2 + C ? % 

lh=\\- —r^-cosf, Ik 


l+4e 2 


l5e 2 + C . 

V^ sm ^ 


(A30) 


(A31) 


for -5e 2 < C < 0 (polar precession). In either case, the third com¬ 
ponent of / follows from the normalisation condition |/| = 1 . 

The instantaneous precession rate |d//df| = |0x/| varies along 
the precession paths. It is minimal at the largest value of \l • k\ along 
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the precession path and drops to zero for l = k (the unstable ori¬ 
entations). The instantaneous precession rate becomes maximal at 
I k- 0 (at zero twist 0), when 

|d//df| max = ^ , -t( 1 +4e 2 )singcosg. (A32) 

4(1 +q) z r l 

Thus, the maximum precession rate is much larger for highly ec¬ 
centric than for circular binaries. The largest variation of the pre¬ 
cession rate occurs for precession paths close to the separatrices. 



